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A  2D  multi-physic  in-house-model  has  been  developed  to  analyse  the  performances  of  Solid  Oxide  Elec¬ 
trolysis  Cell  (SOEC)  stack.  This  model  encompasses  a  combined  electrochemical  and  thermal  description 
of  the  electrolyser.  An  analytical  solution  for  multi-species  diffusion  across  the  porous  cathode  has  been 
implemented  in  the  model.  This  numerical  tool  is  useful  to  provide  all  the  important  parameters  of  the 
stack  operation:  distribution  of  temperature,  heat  fluxes,  local  current  densities,  gas  concentrations  and 
overpotentials. 

Simulations  show  that  thermal  equilibrium  of  the  stack  is  strongly  dependant  on  radiative  heat  losses, 
whereas  convective  heat  transfers  are  limited.  In  the  exothermic  operation  mode,  the  cell  warming 
depends  on  the  radiation  efficiency. 

A  parametric  study  has  been  carried  out  to  analyse  SOEC  irreversible  losses.  It  is  found  that  the  anode 
activation  overpotential  is  significant,  whereas  the  polarisation  due  to  the  cathode  activation  remains 
much  more  limited.  Anode  concentration  overpotentials  are  found  to  be  insignificant  whatever  the  oper¬ 
ating  condition.  Cathode  concentration  overpotentials  are  found  to  be  moderate  in  the  case  of  Electrolyte 
Supported  Cell  (ESC),  while  they  can  increase  drastically  with  current  density  in  the  case  of  Cathode  Sup¬ 
ported  Cell  (CSC).  Addition  of  diluent  gases  into  the  H2,  H20  mixture  is  found  to  increase  the  concentration 
overpotentials. 

©  2010  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

In  the  context  of  fossil  fuel  rarefaction,  hydrogen  is  consid¬ 
ered  as  a  promising  energy  carrier  liable  to  complement  or  even 
replace  fossil  hydrocarbon  fuel  [1].  When  driven  by  renewable  or 
nuclear  energy,  water  electrolysis  allows  producing  “clean”  hydro¬ 
gen  without  greenhouse  gas  emission.  The  energy  required  for 
low  temperature  water  electrolysis,  either  with  acidic  or  alka¬ 
line  electrolyte,  has  to  be  entirely  provided  as  electricity.  For  high 
temperature  water  electrolysis,  namely  High  Temperature  Steam 
Electrolysis  (HTSE),  the  energy  required  is  lower  because  the  water 
molecules  are  already  gaseous.  Moreover  part  of  this  energy  can 
be  brought  as  heat  [2].  As  a  result,  when  heat  sources  are  avail¬ 
able  HTSE  stands  as  a  highly  relevant  hydrogen  production  mode 
[3-5].  Furthermore,  the  high  temperature  allows  avoiding  the  use 
of  expensive  catalysts  for  the  electrochemical  reactions. 

The  steam  electrolysis  reaction  (1)  is  the  reverse  operation  of 
Solid  Oxide  Fuel  Cells  (SOFC).  At  the  cathode  side,  water  is  reduced 
according  to  the  electrochemical  reaction  (la)  written  in  the  case 
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of  an  anion  conducting  electrolyte,  whereas  oxygen  is  produced  at 
the  anode  side  through  reaction  (lb): 

H20  -*■  l/202+H2  AH°(800°C)  =  +242 KJ mol-1  (1) 

H20  +  2e"^  H2  +  02“  (la) 

02~  ->•  1  /202  +  2e~  (lb) 

The  feasibility  of  the  electrolysis  operation  has  been  demon¬ 
strated  in  the  80s  by  Dornier  using  tubular  cells  [6].  More  recently, 
planar  SOFC  cells,  SRUs  and  short  stacks  operated  in  the  electrol¬ 
yser  mode  have  given  promising  hydrogen  production  rates  [7,8]. 
In  these  studies,  Solid  Oxide  Electrolysis  Cells  (SOEC)  are  consti¬ 
tuted  by  the  same  ceramic  assembly  as  SOFC:  the  dense  electrolyte 
is  classically  made  of  Yttria  Stabilised  Zirconia  (YSZ),  porous  Lan¬ 
thanum  doped  Strontium  Manganite  (LSM)  is  used  for  the  anode 
and  a  porous  cermet  composed  of  YSZ  and  Nickel  (Ni-YSZ)  is  used 
for  the  cathode.  In  planar  configuration,  cells  are  stacked  into  a 
high  packing  density  and  are  separated  by  stainless  steel  plates, 
ensuring  the  electrical  inter-connection  and  gas  distribution  on  the 
electrodes.  However,  if  high  performance  levels  have  been  obtained 
by  operating  SOFC  in  the  electrolyser  mode,  some  limitations  have 
also  been  highlighted  regarding  SOEC  durability  either  at  the  single 
cell  level  [8]  or  when  integrated  into  a  stack  [9,10]. 
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List  of  symbols 

a  manifold  width  (cm) 

b  manifold  height  (cm) 

L  manifold  length  (cm) 

Clp  heat  capacity  (J  mol-1  K-1 ) 

Defl  effective  molecular  diffusion  coefficient  (cm2  s-1 ) 

Dfjc  effective  Knudsen  diffusion  coefficient  (cm2  s-1 ) 

effective  diffusion  coefficient  (cm2  s_1 ) 

Ei= o  Open  Circuit  Voltage  (mV) 

F  Faraday  constant  (=96485  C  mol-1 ) 

h  heat  transfer  coefficient  (W  nrr2  K-1 ) 

i  local  current  density  ( mA  cm-2 ) 

z'o  exchange  current  density  ( mA  cm-2 ) 

rii  molar  flow  rate  (mol  s-1 ) 

Nj  molar  flux  (mol  s-1  cm-2 ) 

r  mean  pore  radius  (p,m) 

R  gas  constant  (=8.314J  mol-1  K-1 )  or  H2/H20  ratio 

Rohm  ohmic  resistance  ( £2  cm2 ) 

Rc  contact  resistance  ( £2  cm2 ) 

Rlsm  electrical  resistance  of  anode  ( £2  cm2 ) 

RNi-8Ysz  electrical  resistance  of  cathode  ( £2  cm2 ) 

T  temperature  (I<) 

U  cell  voltage  (mV) 

yi  molar  fraction 

Sa  anode  thickness  (cm) 

8C  cathode  thickness  (cm) 

8e  electrolyte  thickness  (cm) 

A Hr  reaction  enthalpy  (J  mol-1 ) 

X  thermal  conductivity  (W  nrr 1  K-1 ) 

r]  dynamic  viscosity  (kg  nrr1  s-1 ) 

77act  anode  or  cathode  activation  overpotentials  (mV) 

Tjconc  anode  or  cathode  concentration  overpotentials  (mV) 

s  porosity  or  emissivity 

r  tortuosity  factor 

p  electrical  conductivity  ( cm-1 ) 

( 7  Stefan-Boltzmann’s  constant  (=5.67  x  1 0-8  W nr2  K-4) 

§  cathodic  gas  dilution  rate 

Computing  simulation  appears  to  be  one  of  the  most  efficient 
approaches  to  analyse  the  coupled  mechanisms  of  SOEC  opera¬ 
tion.  Indeed,  the  cell  response  depends  on  multi-physic  phenomena 
which  are  entangled  such  as  the  electrochemical  reactions,  mass 
and  energy  transport  in  the  electrolyser.  Some  models  have  been 
recently  developed  to  describe  the  SOEC  behaviour.  Ni  et  al.  [  1 1 ,1 2  ] 
have  proposed  isothermal  models  taking  into  account  the  coupled 
mass  transport  and  electrochemical  reactions.  In  reference  [12], 
diffusion  of  H20  and  EI2  species  at  the  cathode  side  are  described 
through  the  Dusty  Gas  Model  (DGM).  Udagawa  et  al.  [13,14]  have 
developed  a  one  dimensional  dynamic  model  including  a  thermal 
analysis.  For  design  optimisation,  commercial  Computational  Fluid 
Dynamic  (CFD)  software  has  been  used  by  Herring  et  al.  [  1 5 ]  to  pro¬ 
pose  a  three  dimensional  SOEC  modelling.  Freshly,  Gondrin  et  al. 
[16]  have  improved  the  two  dimensional  CFD  description  by  imple¬ 
menting  the  DGM  approach  to  calculate  H2  and  H20  distribution 
through  porous  cathode  as  a  function  of  microstructure. 

In  the  present  work,  an  in-house  model  has  been  developed  to 
describe  the  thermo-electrochemical  behaviour  of  a  planar  SOEC 
stack.  A  two  dimensional  approach  has  been  adopted  in  order  to 
couple  locally  fluidic,  thermal  and  electrochemical  steps.  This  SOEC 
model  is  derived  from  a  first  one  previously  developed  for  SOFC 
operation  simulation  [17,18].  The  mass  transfer  into  the  porous 
electrodes  has  been  described  in  the  frame  of  the  DGM  model.  At  the 
cathode  side,  the  effect  of  gas  dilution  has  been  taken  into  account 
by  adding  nitrogen  into  the  H20/H2  mixture.  Mass  transport  and 
kinetics  of  electrochemical  reactions  have  been  coupled  to  a  ther- 


Fig.  1.  Schematic  representation  of  the  simulated  SRU  (cross-section  view)  and 
coordinate  system  used  in  the  model. 

mal  analysis.  A  special  attention  has  been  paid  on  radiative  heat 
losses  that  are  usually  neglected  despite  their  significant  influence. 
Mechanisms  of  steam  electrolysis  and  cell  operation  have  been 
analysed.  The  effects  of  H2/H20  ratio  introduced  at  the  cathode 
inlet  and  the  impact  of  gas  dilution  on  the  electrolyser  efficiency 
have  been  investigated. 

2.  Model  description 

2.1.  Introduction:  geometry  and  materials 

A  schematic  cross-section  of  the  simulated  Single  Repeat  Unit 
(SRU)  is  illustrated  in  Fig.  1 .  The  cell  presents  a  square  shape  and  the 
electrodes  are  centred  on  the  electrolyte.  Most  classically  materials 
used  for  SOEC  have  been  considered  for  the  simulations :  LSM  for  the 
anode,  Zr02  stabilised  with  8  mol%  Y203  (8YSZ)  for  the  electrolyte 
and  Ni-8YSZ  cermet  for  the  cathode. 

Interconnects  are  assumed  to  be  made  in  ferritic  stainless  steel 
plates.  A  glass  seal  is  positioned  into  a  groove,  which  is  machined 
in  the  metallic  plate  in  order  to  prevent  the  gas  leakage  (Fig.  1 ). 

Air  and  steam  are  introduced  respectively  at  the  anode  and  cath¬ 
ode  sides  by  two  gas  manifolds  designed  into  the  interconnects  and 
located  at  the  electrodes  gas  inlet  (Fig.  1).  These  manifolds  ensure 
gas  feeding  into  a  pattern  of  triangular  channels.  As  illustrated  in 
Fig.  2,  these  channels  distribute  the  gases  on  each  electrode  in  a 
counter-flow  configuration.  The  produced  hydrogen  and  oxygen 
are  finally  exhausted  from  the  SRU  and  collected  by  the  two  out¬ 
let  gas  manifolds,  which  are  designed  at  the  electrodes  gas  outlets 
(Fig.  1).  The  complete  gas  flow  scheme  (taken  here  for  air)  is  rep¬ 
resented  in  Fig.  3.  Some  typical  dimensions  of  the  SRU  are  listed  in 
Table  1. 

The  coordinate  system  used  for  the  modelling  is  shown  in 
Figs.  1  and  3.  A  two  dimensional  approach  has  been  adopted.  The 
modelled  surface  is  contained  in  the  (x,y)  plane  (perpendicular  to 


Fig.  2.  Configurations  for  the  oxygen  and  hydrogen  flows  along  the  electrodes 
(counter-flow  configuration).  The  triangular  shape  of  the  anode  and  cathode  gas 
channels  is  detailed. 
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the  z  direction)  and  corresponds  to  an  elementary  slice  taken  in 
the  middle  of  the  SRU  along  a  gas  channel.  Therefore,  this  approach 
neglects  edge  effects  on  the  z  direction. 

The  model  has  been  divided  into  two  connected  modules: 
the  electrochemical  one  and  the  thermal  one.  The  electrochem¬ 
ical  part  of  the  model  has  been  implemented  in  the  software 
tool  Matlab  [19].  It  allows  determining  the  local  current  density, 
the  distribution  of  gas  species  and  the  thermal  sources  related 
to  the  electrochemical  reactions.  The  thermal  module  has  been 
implemented  in  a  finite  element  solver,  Cast3M  [20].  It  allows 
determining  the  thermal  fluxes  and  temperature  field  in  both  the 
solid  and  gas  phases  of  the  SRU.  The  two  modules  have  been  cou¬ 
pled  through: 

(i)  the  temperature  field  required  to  determine  the  kinetic  con¬ 
stants  of  electro-chemical  reactions, 

(ii)  the  heat  source  terms  necessary  to  simulate  the  temperature 
field, 

(iii)  the  gas  composition  variation  along  the  channels  necessary  to 
calculate  the  fluids  thermal  properties. 


the  DGM  model,  by  combining  the  Stephan-Maxwell  and  Knudsen 
diffusion: 


Ni  A  yA-yiNj  PT  /  dy,  \ 

Deff  +  2^  Deff  RT  \dy  J 

Uk  Ui 


where  yx  is  the  molar  fraction  of  specie  i  and  Pj  the  total  pressure 
(=1  atm).  The  effective  Knudsen  and  binary  diffusion  coefficients, 
D^k  and  Df-J,  have  been  determined  as  a  function  of  the  electrodes 
microstructure  parameters  (i.e.  mean  pore  radius  r,  tortuosity  fac¬ 
tor  r  and  porosity  s).  In  porous  media  such  as  SOFC  or  SOEC 
electrodes,  the  following  expressions  are  usually  employed  [24,25 ] : 


4!  -  SrDUi  and  DU  =  7D>.* 


(4) 


The  Knudsen  diffusion  coefficient  can  be  assessed  according  to  the 
kinetic  theory  of  gases: 


(5) 


where  Mz  denotes  the  molecular  weight  for  the  gas  t.  The  binary 
diffusion  coefficients  have  been  expressed  according  to  Todd  and 
Young  [26]: 


DiJ  = 


0.00143T1-75 


PrMlfivf3  +  V‘/3) 


V3\ 


With  My  =  2/(  T  +  jl 


(6) 


where  corresponds  to  the  Fuller  diffusion  volume. 

At  the  anode  side,  the  nitrogen  flux  is  equal  to  0  (NN2  =  0)  since 
this  component  is  neither  consumed  nor  produced.  The  oxygen 
release  rate  at  the  anode/electrolyte  interface  is  given  by  the  local 
current  density  i  according  to  the  Faraday’s  law  (note  that  i  is  con¬ 
sidered  as  its  absolute  value): 


2.2.  Description  of  the  electrochemical  module 
2.2.1.  Mass  transfer 

The  pressure  drop  along  the  gas  channels  has  been  assumed  to 
be  negligible:  atmospheric  condition  has  been  considered  for  the 
gas-phase  pressure  (PT  =  1  atm). 


No2  =  jp  (for  y  =  electrolyte  thickness)  (7) 

Considering  these  conditions,  the  method  used  in  [21  ]  for  inte¬ 
grating  (3)  in  SOFC  mode  can  be  also  applied  for  SOEC  mode 
and  leads  to  the  expression  of  the  oxygen  molar  fraction  at  the 
anode/electrolyte  interface  yg£: 


2.2.1. 1.  Gas  flow  along  the  anodic  and  cathodic  channels.  The  varia¬ 
tion  of  molar  fractions  along  the  gas  channels  (i.e.  in  the  x  direction) 
is  caused  by  molar  fluxes  through  the  porous  electrodes  (i.e.  in  they 
direction).  This  variation  is  calculated  by  using  local  mass  balances 
performed  for  each  species  along  the  channels: 


where  n2  denotes  the  molar  flow  rate  in  the  channels  for  the  specie 
i,  Nj,  the  flux  through  the  electrode,  e,  the  width  of  the  gas  channel 
and,  yi,  a  stoichiometric  coefficient  (equal  to  -1  for  the  reactant 
species  and  +1  for  the  product  species). 

2.2.1. 2.  Mass  transport  in  the  porous  electrodes.  In  porous  elec¬ 
trodes,  the  viscous  flow,  which  is  driven  by  a  pressure  gradient,  is 
generally  considered  as  negligible  compared  to  the  diffusive  flow 
[21-23].  Therefore,  the  mass  transport  through  the  porous  elec¬ 
trodes  (i.e.  in  the  y  direction)  has  been  described  in  the  frame  of 


Table  1 

Typical  dimensions  of  the  simulated  SRU. 


Electrodes 

Manifold 

Manifold 

Manifold 

Interconnect 

length,  l 

length,  L 

height,  h 

wide,  b 

plate  thickness, 

(mm) 

(mm) 

(mm) 

(mm) 

w  (mm) 

88 

91 

3 

3 

10 

y'Sl  =  1  +  (yCo',Ml  ~  1  )exP  (  — RT!f  )  (for  PT  =  1  atm  and  y 

2  2  V  4 FDf2PT ) 

=  electrolyte  thickness)  (8) 

where  y™™1  is  the  oxygen  concentration  in  the  channel  and  8a  the 
anode  thickness. 

At  the  cathode  side,  the  nitrogen  molar  flux  is  also  equal  to  0 
(NN2  =  0).  The  fluxes  related  to  the  water  consumption  and  hydro¬ 
gen  production  at  the  cathode/electrolyte  interface  are  also  given 
by  the  Faraday’s  law: 

Nh2o  =  -Nh2  =  2p  (for  pT  =  1  atm  and  y  =  0)  (9) 

A  methodology  detailed  in  Appendix  A  is  proposed  here  to  inte¬ 
grate  the  system  of  equations  (3)  in  the  case  of  a  gas  mixture 
composed  of  hydrogen,  steam  and  nitrogen  (used  as  a  dilution  gas). 
It  allows  determining  the  molar  fractions  at  the  cathode/electrolyte 
interface: 

yfr*  =  y™"al  x  exP  (-fci  sc)  (i  oa) 

<  =  VhT'  +  Me  +  |  x  (y“"“(  -  exp  (-Me)  )  (10b) 

y^o  =yH2o'  +  Mc  +  |  x  (yj?™'  -yj^'exp  (-Me))  (10c) 
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Microstructure  and  electrochemical  parameters  of  cell  layers. 
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Porosity,  e 

Tortuosity  factor,  r 

Mean  pore  radius,  r  (p,m) 

Microstructural  parameters  of  both  electrodes  (anode  and  cathode) 

0.4  4 

Global  contact  resistances  between  electrodes/interconnects :  Rc  =  0.1  £2  cm2  [27] 
Electrodes  electronic  conductivities :  Pm-ysz  =  800  cm-1  [33,34];  Plsm  =  72  cm- 

i 

1  [35] 

Pre-exponential  factors  (mAcm-2) 

Molar  fraction  exponent 

Activation  energy  (KJ  mol-1 ) 

Exchange  current  densities 

Anode  (LSM) 

Cathode  (Ni-8YSZ) 

B  =  2.05x  108 

A  =  1.26  x  1010 

p  =  1/4  [36,37] 
m  =  n  =  l  [36,38] 

120  [36] 

120  [36] 

where  ki= i_5  are  constants  given  in  Appendix  A,  that  depends  on 
diffusion  coefficients,  local  temperature  and  current  density. 

2.2.2.  Cell  voltage 

Electrodes  have  been  assumed  to  be  good  electrical  conductors 
(Table  2)  so  that  the  electrical  potential  can  be  considered  constant 
along  the  cell.  Therefore,  in  the  present  model,  the  cell  voltage  is 
fixed  and  the  local  current  density  is  calculated  in  order  to  verify 
the  following  equation  along  the  cell  length: 

ii  r  y.anode  ,  „ cathode  ,  „anode  ,  M cathode  ,  n  :  nu 

ucell  =  hi=0  +  Bconc  +  Bconc  +  Bact  +  Vact  +  Kohml  l 1  1 ) 


be  expressed  as  a  function  of  the  current  density  from  the 
Butler-Volmer  equation: 

(15a) 

{.»p  -,-p  (_2(1  -  ’Jpsri) }  (i5b) 

where  io<anode  and  io, cathode  represent  respectively  the  anode  and 
cathode  exchange  current  densities.  With  the  symmetrical  factor 
a  taken  equal  to  0.5,  the  activation  overpotentials  can  be  derived 
from  Eqs.  (15a)  and  (15b)  [30,31]: 


with  Ei= o  corresponds  to  the  Open  Circuit  Voltage  (OCV)  expressed 
from  the  Nernst’s  equation: 


nT  Vo 

Ei=0  =  E°(T)+^  In 


,i=0  (yTPBa,i=0\0  5 


y02 
VTPBCJ= 

yH  20 


(12) 


where  E°  is  the  standard  potential  and  yj  the  partial  pressures  at 
OCV  (i  =  0)  or  under  current  (i  =f=  0).  The  superscripts  TPBa  and  TPBC 
denote  the  Triple  Phase  Boundaries  at  anode  and  cathode  sides. 

The  terms  rjCOnc  and  pact  in  Eq.  (11)  are  expressed  considering 
their  absolute  values  and  correspond  respectively  to  the  concentra¬ 
tion  and  activation  overpotentials  (at  the  anode  or  cathode  side). 
The  term  R0hmi  denotes  the  pure  ohmic  losses  including  (i)  the  ionic 
resistance  of  the  electrolyte,  Re,  (ii)  electronic  resistances  of  both 
electrodes,  Rlsm  and  Rm-sysz*  and  (iii)  the  global  contact  resistances 
between  the  electrodes  and  the  interconnects,  Rc: 


Rohm  —  R-LSM  +  RNi-8YSZ  +  Re  +  Rc 


(13a) 


with  (S  =  1  cm2 )  RLSm  = 


8q 

Plsm  ’ 


*Ni-8YSZ  = 


PNi-8YSZ 


Re  = 


P8YSZ 


(13b) 


where  p  is  the  electrical  conductivity  of  each  considered  material 
and  8  the  thickness  of  the  anode,  electrolyte  or  cathode  layer.  The 
contact  resistance  Rc  is  due  to  the  current  collection  between  elec¬ 
trodes  and  interconnects.  The  value  of  this  parameter  has  been 
taken  in  the  model  to  0.1  £2  cm2  [27]  (Table  2).  In  planar  config¬ 
uration,  electronic  resistance  of  the  electrodes  is  low  compared  to 
ohmic  losses  due  to  the  electrolyte  and  can  be  regarded  as  temper¬ 
ature  independent  within  the  operating  conditions  [28].  Inversely, 
the  ionic  conductivity  of  the  8YSZ  electrolyte  has  to  be  expressed 
as  function  of  the  operating  temperature  and  is  taken  as  [29]: 


pe  (£2  1  cm  1 )  =  466  x  exp 


(14) 


2.2.2 A.  Activation  overpotentials.  The  activation  overpotentials 
represent  the  voltage  loss  induced  by  the  electrochemical  reac¬ 
tions  occurring  at  the  electrodes/electrolyte  interfaces.  They  can 


RT  •  u-l 
+  —  sinh 
F 


L  | 

^  ^  x  *0,  cathode  J 


(16) 


The  exchange  current  densities  depend  on  (i)  the  temperature 
according  to  an  Arrhenius  law  and  (ii)  to  the  reactant/product  con¬ 
centrations  in  the  neighbourhood  of  the  active  sites  (i.e.  the  TPBs) 
[32]: 

Chode  =A  X  (y£f  (<  =  0))m  X  (y™c(i  =  0))"  X  e(-^tW7KT)  (17a) 
/™ode  =  B  x  (y™a(i  =  0))p  x  e(-Code/RT)  (17b) 

The  values  chosen  for  the  activation  energies  Ea  and  the  molar 
fraction  exponents  are  reported  in  Table  2.  These  values  have  been 
determined  on  typical  Ni-YSZ//YSZ//LSM  cells  operated  in  SOFC 
mode.  The  pre-exponential  factors  A  and  B  have  been  scaled  in  order 
to  obtain  i0  values  recommended  in  literature  [31,39].  It  is  have 
assumed  in  this  work  that  these  exchange  current  density  expres¬ 
sions  (17a)  and  (17b)  were  reasonably  applicable  to  SOEC  operation 
since  temperature  and  reactant/product  species  are  similar  in  both 
SOFC/SOEC  operation  modes. 


2.2.22.  Concentration  overpotentials.  Potential  losses  can  arise  in 
operation  because  of  the  electrode  inability  to  maintain  the  fluid 
initial  composition  in  the  gas-phase  surrounding  the  TPBs.  These 
concentration  overpotentials  are  expressed  from  the  Nernst’s  equa¬ 
tion.  They  can  be  caused  by  insufficient  diffusion  through  the 
electrodes  or  insufficient  steam  flow  rate  introduced  at  the  SOEC 
inlet: 


y.cathode 

'leone 


RT  ,  yj&(l  =  0)x3^(f#0) 
2 F  \>V0)xy^(i  =  0) 


anode  .  RT  0) 


Vconc 


—  4E  ny™a(i  =  0) 


(18a) 

(18b) 


It  can  be  noticed  that  the  molar  fraction  at  the  TPBa  and  TPBC 
have  been  reduced  in  the  present  model  at  the  anode/electrolyte 
and  cathode/electrolyte  interfaces  respectively.  It  is  worth  noting 
that  this  assumption  is  well  verified  for  sufficiently  thick  electrodes. 
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Fig.  4.  Boundary  conditions  assumed  for  the  thermal  simulations. 


Indeed,  several  studies  [21,40,41]  have  shown  that  the  reaction 
zone  spreads  on  a  limited  region  from  the  electrolyte/electrode  into 
the  electrode  (~10  |jim). 

2.3.  Description  of  the  thermal  module 

2.3 A.  Boundary  conditions 

Adiabatic  conditions  have  been  assumed  for  the  surfaces  con¬ 
necting  two  adjacent  SRUs  (i.e.  the  bottom  and  top  free  surfaces  of 
SRU  as  shown  in  Fig.  4).  This  assumption  is  rigorously  right  only  for 
the  SRUs  located  in  the  central  region  of  the  stack  where  heat  flux 
in  the  stacking  direction  can  be  neglected. 

The  temperature  of  gases  introduced  into  the  SRU  (i.e.  at  the 
two  gas  manifold  inlets:  see  Fig.  3)  have  been  taken  equal  to  800  °C. 
The  temperature  of  the  insulating  envelope  surrounding  the  stack 
is  also  supposed  to  be  maintained  at  800 °C.  This  last  assumption 
corresponds  to  either  a  stack  tested  in  a  furnace,  or  a  system  in 
which  the  thermal  fluxes  exchanged  with  the  SRU  are  managed. 

2.3.2.  Solid  temperatures:  energy  balance  for  solid  phases  into  the 
SRU 

The  governing  equation  to  calculate  the  temperature  field  in 
the  solids  has  been  expressed  considering  conduction,  convection 
and  radiation.  For  a  solid  volume,  dV ,  in  contact  with  the  fluid  on  a 
surface,  dS,  energy  balance  has  been  written  as  follows: 

A  {div(grad(T))}  dV  +  dQ  =  hgasdS  {Ts-Tg}+  d<prad 

t  t  t  t  (19) 

(conduction)  +  (heat  sources)  =  (convection)  +  (radiation) 

The  first  term  of  this  equation,  X  { div(grad(T ))}  dV,  is  related  to  the 
heat  transport  by  conduction  in  the  solid  phases.  The  term  X  rep¬ 
resents  the  thermal  conductivity  of  the  solid  volume  dV.  It  is  worth 
noting  that  the  porous  anode  and  cathode  have  been  modelled  as 
homogenous  media  in  which  only  conduction  is  taken  into  account. 
Indeed,  the  Peclet  number  Pe  calculated  within  both  electrodes 
remains  much  lower  than  unity,  meaning  that  the  heat  transport 
by  convection  is  negligible  in  comparison  to  conduction  [42].  The 
effective  thermal  conductivity  of  the  porous  electrodes  can  then  be 
determined  considering  a  mixture  law  between  the  conductivities 
of  solid  and  fluid  regions  [43]. 

2.3.3.  Sources  terms 

The  second  term  of  Eq.  (19),  dQ,  corresponds  to  the  heat  sources 
related  to  the  electrochemical  reaction: 

=  {+w AHr  ~  iMU«*}ds  (20) 

where  A Hr  denotes  the  enthalpy  of  reaction  for  the  steam  disso¬ 
ciation  into  hydrogen  and  oxygen  (see  Eq.  (1)).  The  heat  sources 
include  several  contributions  corresponding  to  the  ohmic  dissipa¬ 
tions  due  to  contact,  electrodes,  electrolyte  resistances  and  heat 


absorption  due  to  the  electrochemical  reaction.  Location  within  the 
SRU  of  these  thermal  sources  is  detailed  in  reference  [18]. 

2.3.4.  Gas  temperatures  for  convective  heat  exchanges 

Computation  of  the  convective  term  in  Eq.  (19)  requires  the 
knowledge  of  the  fluids  temperatures  Tg.  These  fluids  tempera¬ 
tures  must  be  determined  all  along  the  electrode  lengths  in  the 
triangular-shaped  channels  (see  Section  2.3.4.1)  and  (ii)  into  the 
manifolds  at  the  inlet  and  outlet  part  of  the  electrodes  (see  Section 
2.3.42). 

2.3.4A.  Energy  balance  for  fluids  along  the  gas  channels  (in  the 
x  direction).  Because  of  the  electrochemical  reactions,  the  gas 
composition  is  modified  along  the  triangular-channels  between 
interconnect  and  electrodes.  Therefore,  the  gas  properties  (i.e. 
heat  capacity  Cp,  thermal  conductivity  X,  and  dynamic  viscosity  rj) 
depend  on  the  position  x  along  the  channels.  In  this  work,  they 
have  been  calculated  with  the  polynomial  curves  established  by 
Todd  and  Young  [26]. 

Gas  inlet  flow  rates  have  been  chosen  in  such  way  that  a  lami¬ 
nar  flow  takes  place  inside  the  channels  (in  the  studied  conditions, 
it  has  been  found  that  the  Reynolds  number  Re  remains  much 
lower  than  1000  everywhere  in  the  ducts  of  the  simulated  SRU). 
Furthermore,  according  to  the  results  of  Damn  and  Fedorov  [44] 
and  Sanchez  et  al.  [45],  the  radiative  heat  exchange  between  gases 
and  walls  can  be  neglected.  Finally,  along  the  simulated  triangular¬ 
shaped  channels,  the  energy  balances  have  been  written  as  follows : 

^  (n/ (x)CTp (x)T^r(x))  ^  _  h  ,vVJC  fT  M  T 

/  - - dx  —  hsteam{x)dS  |Ts(x)  —  Tg(x)| 

i 

(i  =  H2,  H20,  N2)  (21a) 

v— '  d  (niWCp(x)Tg(x)) 

2  y— - fa - ~dx  =  KirWdS  {Ts(x)  -  Tg(x)} 

i 

(i  =  02,N2)  (21b) 

The  term  Tg  denotes  the  temperature  of  fluids  and  Ts  denotes 
the  wall  temperature  of  the  solid  phase.  The  exchange  surface 
dS  corresponds  to  the  contact  surface  between  gas  and  solid  for 
the  considered  triangular-shaped  channel.  The  heat  transfer  coef¬ 
ficients  hsteam  and  hair  have  been  determined  from  the  Nusselt 
number  Nu: 

hsteam  or  air (x)  =  i  =  02//N2  Or  i  =  H20//H2//N2 

(22) 

where  DH  is  the  hydraulic  diameter  of  the  cathode  or  anode  chan¬ 
nel.  The  determination  of  Nu  number  is  potentially  complex  since 
it  could  be  affected  by  many  factors  such  as  thermal  gradient  or 
heat  flux  heterogeneity.  However,  the  thermal  gradient  remains 
limited  in  the  present  study  (see  Section  3.1.2).  Furthermore,  the 
Nusselt  number  tends  to  an  asymptotic  value,  provided  that  a  fully 
developed  laminar  flow  is  achieved  in  the  channels.  This  condition 
is  largely  fulfilled  in  our  configuration.  As  a  consequence,  the  Nus¬ 
selt  number  could  be  taken  to  Nu  =  3.68  (i.e.  the  asymptotic  value 
corresponding  to  a  flow  confined  in  channel)  [46,47]. 

It  is  worth  noting  that  fluids  heat  conductivities  X[  depend  on 
the  gas  composition  along  the  cell  (Fig.  5).  As  a  consequence,  the 
heat  transfer  coefficient,  which  is  linearly  dependant  of  Xt  accord¬ 
ing  to  Eq.  (22),  will  follow  the  same  evolution.  In  electrolysis  mode, 
steam  is  transformed  into  hydrogen  all  along  the  cathode.  Since 
this  gas  component  presents  a  much  higher  thermal  conductivity 
than  other  species,  the  heat  exchange  coefficient  hsteam(x )  can  be 
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Fig.  5.  Evolution  of  the  fluid  thermal  conductivity  and  heat  transfer  coefficient  plot¬ 
ted  along  the  cathode  channel.  Nu  number  has  been  assumed  constant  along  the 
channel  (Uceii  =  1.5  V,  steam  conversion  rate  =  96.2%,  condition  I  of  SRU  gas  feeding: 
see  Table  4,  active  cell  area  =  77.44  cm2 ). 


strongly  increased  from  the  electrode  inlet  up  to  the  outlet  (Fig.  5). 
From  this  result,  it  appears  that  modification  of  the  heat  transfer 
coefficient  along  the  cell  has  to  be  taken  into  account  for  SOEC 
modelling. 

23.42.  Fluids  temperature  into  the  manifolds.  Convective  heat 
exchanges  have  been  taken  into  account  for  gas  flowing  into  the 
manifolds: 


•  The  first  one,  d(p\ ,  corresponds  to  the  solid  surface-to-surface  heat 
exchange  between  electrodes  (anode  or  cathode)  and  intercon¬ 
nects  plates.  It  has  been  approximated  in  the  model  according  to 
the  general  expression  between  two  infinite  parallel  planes: 


dip\  =  • 


&£  anode  I  cathode  ^interconnect 


1  ^1  £  anode  I  cathode')  (1  £  interconnect ) 


74 

anode/cathode 


T4  ) 

interconnect  / 


.  dS 


(24) 


where  s  denotes  the  emissivity  of  materials  (for  anode,  cathode 
and  interconnects)  and  o  the  Stefan-Boltzmann’s  constant. 

•  The  second  contribution  to  diprad  corresponds  to  the  radiative 
heat  losses,  dcp2,  flowing  from  the  SRU  to  the  stack  insulating 
envelope: 

d<P2  =  {^interconnect  (js  ~  ^insulating)  }  (25) 

where  Ts  denotes  the  surface  temperature  taken  on  the  free  edge 
of  the  SRU.  The  term  Tinsuiating  corresponds  to  the  stack  envelope 
temperature  taken  here  at  800  °C  (Fig.  4).  The  thermal  parameters 
for  materials  considered  in  the  model  are  provided  in  Table  3. 


The  model  reliability  can  be  considered  as  good  since  it  is  based 
on  usual  physics  for  the  studied  SOEC  application.  Furthermore,  the 
model  derived  from  a  SOFC  one  [17,18]  which  has  been  already  val¬ 
idated.  It  can  be  noticed  that  experimental  validation  of  the  present 
model  is  on  going  and  will  be  published  later. 


•  At  the  electrodes  inlets,  the  calculation  of  steam  (or  air)  temper¬ 
ature  evolution  in  the  z  direction  is  provided  in  Appendix  B.  The 
convective  heat  fluxes  introduced  in  Eq.  (19)  have  been  deter¬ 
mined  considering  the  gas  temperature  taken  in  the  middle  of 
the  SRU  ( i.e .  for  z  =  L/2:  see  model  assumptions  in  Section  2.1): 

r  i  'j  ( h2(a+b)L/nt0tCp ) 

Tg  =  (Tiniet-Ts)x 

+  TS  (for  z=l-  and  Tinlet  =  800 °c)  (23) 

where  Tg  is  the  steam  (or  air)  temperature  into  the  manifold.  Tin/et 
and  ntot  are  respectively  the  fluid  temperature  and  molar  flow 
rate  introduced  in  the  SRU. 

•  At  the  electrodes  outlet,  the  gas  temperature  into  the  manifold 
has  been  assumed  to  be  the  same  that  the  one  taken  at  the  outlet 
of  the  triangular-shaped  channels  (i.e.  for  x  =  l  and  z  =  L/2). 

2.3.5.  Radiative  heat  exchanges 

The  last  term  of  Eq.  (19),  d(prad ,  is  related  to  the  heat  transfer  by 
radiation  and  has  been  divided  into  two  contributions: 


3.  Simulation  results  and  analysis 

3.1.  Simulations  of  a  case  taken  as  reference:  mechanisms 
involved  in  the  steam  electrolysis 

3.1.1.  Simulation  conditions  for  the  case  taken  as  reference 

A  reference  case  has  been  defined  in  order  to  analyse  the  mech¬ 
anisms  involved  in  the  steam  electrolysis.  In  this  reference  case,  an 
Electrolyte  Supported  Cell  (ESC)  with  anode,  electrolyte  and  cath¬ 
ode  thicknesses  respectively  fixed  to  50, 90  and  50  pan  is  integrated 
in  the  SRU.  The  cell  active  area  is  taken  to  77.44  cm2. 

Simulations  have  been  carried  out  for  the  first  condition  of  gas 
feedings  detailed  in  Table  4.  These  conditions  have  been  chosen 
to  achieve  a  high  steam  conversion  (SC)  rate  at  Uceu  =  1.5  V  in  such 
way  that  the  mass  and  energy  transports  phenomena  associated  to 
the  overpotentials  are  exacerbated.  It  can  be  noticed  that  the  same 
steam  and  hydrogen  molar  flow  rate  is  supplied  at  the  cathode  side : 
R  defined  as  the  ratio  of  nH2  over  nH2o,R  =  nH2  /nH2o,  is  equal  to  the 
unity.  Moreover,  the  dilution  of  hydrogen  and  steam  into  nitrogen 
is  limited:  §  expressed  as  §  =  nN2/(nH2  +  nH2o  +  nN2)  is  taken  to  0.1. 

Air  ratio  X  is  defined  as  the  air  molar  flow  rate  introduced  at  the 
anode  inlet  divided  by  the  cathode  molar  flow  rate  (A  =  nair/nca^igas ). 


Table  3 

Thermal  parameters. 


Cathode 

Electrolyte 

Anode 

Interconnect 

Glass  seal 

Conductivity,  X  (W  irr 1  K_1 ) 

6  [43] 

3.8  [48] 

11  [43] 

25  [49] 

1.0a 

Emissivity,  s 

0.4  [43] 

- 

0.4  [43] 

0.7  [50] 

- 

a  Estimated  with  correlations  provided  in  ref.  [51  ]  considering  a  commercial  glass  seal  (Schott  8422). 


Table  4 

Simulated  molar  flow  rates  (active  cell  area  =  77.44 cm2). 


Gas  feeding  condition 

Cathode  gas  mixture  of  H20,  H2  and  N2 

Anode  gas  (air) 

Hh2o  (xlO  6 mol  s  1  cm-2) 

R  =  nH2/nH2o 

£  =  J1n2  /(nH2  +  Hh20  +  Hn2  ) 

Hair  (  x  1 0-6  mol  S_1  CUT2  ) 

X  -  na\r  In  cathodic  gas 

I  (reference) 

6.9034 

1/1 

1/10 

30.682 

2 

II 

6.9034 

1/9 

1/2 

30.682 

2 

III 

6.9034 

1/1 

1/2 

55.227 

2 
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H2  outlet 
Air  inlet 


Fig.  6.  Evolution  of  the  molar  fraction  along  the  cell  length  in  the  reference 
case  (Uceu  =  1.5  V,  steam  conversion  rate  =  96.3%,  condition  I  of  cell  gas  feeding, 
<5 electrolyte  =  90  (Jim,  8an0de  =  8 cathode  =  50  jjim,  cell  active  area  =  77.44  cm2 ). 


Position  along  the  cell  x  (cm) 

Fig.  8.  Evolution  of  the  overpotential  along  the  cell  length  in  the  reference 
case  (Uce//  =  1.5V,  steam  conversion  rate  =  96.3%,  condition  I  cell  gas  feeding, 
$  electrolyte  =  90  jjim,  8anode  =  8cathode  =  50  jjim,  cell  active  area  =  77.44  cm2 ). 


This  parameter  is  taken  to  2  in  order  to  obtain  a  sufficient  air  flow 
rate,  and  hence,  to  maintain  a  substantial  convective  heat  transfer 
at  the  anode  side. 

3A.2.  Evolution  of  molar  fractions,  overpotentials  and 
temperature  along  the  cell  length 

The  variations  of  the  output  parameters  of  the  model  along  the 
cell  length  are  presented  for  a  steam  conversion  rate  SC  =  96.3%  cor¬ 
responding  to  a  cell  voltage  equal  to  1500  mV.  Fig.  6  illustrates  the 
evolution  of  molar  fractions  calculated  here  at  the  electrochemical 
interfaces.  The  steam  concentration  is  found  to  be  decreased  from 
the  cathode  inlet  up  to  the  cathode  outlet,  since  this  specie  is  con¬ 
sumed  by  the  electrochemical  reaction  ( 1  a).  Accordingly,  hydrogen 
is  continuously  produced  and  its  molar  fraction  is  raised  from  the 
cathode  inlet  up  to  the  outlet.  At  the  anode  side,  the  oxygen  is  pro¬ 
duced  by  reaction  (lb),  and  hence,  its  molar  fraction  increases  from 
the  anode  inlet  up  to  the  anode  outlet. 

Fig.  7  shows  the  evolution  of  the  local  current  density  along 
the  cell  length.  The  current  slightly  decreases  from  the  anode  inlet 
up  to  the  outlet.  This  trend  is  consistent  with  the  steep  increase 
of  the  cathode  concentration  overpotential  pccf^c  calculated  along 
the  cell  length  and  shown  in  Fig.  8.  This  pronounced  increase  can 
be  attributed  to  the  major  water  depletion  associated  to  the  high 
steam  conversion  rate  (SC  =  96.3%). 

It  is  worth  noting  that  the  H20  or  FI2  molar  fraction  gra¬ 
dients  along  the  y-axis  between  the  gas  channel  and  the 
cathode/electrolyte  interface  remains  negligible  all  along  the  cell 


H 2  outlet 
Air  inlet 


length  ( AyH2  0r  h2o  <  0.01).  This  result  means  that  the  gas  trans¬ 
fer,  either  steam  or  hydrogen,  through  the  thin  cathode  is  not  a 
limiting  step  for  the  electrochemical  reaction. 

As  shown  in  Fig.  8,  the  ohmic  voltage  losses  Rohm  x  i  and  the 
anode  activation  overpotentials  paJ^de  are  found  to  follow  the  same 
evolution  than  the  current  density  along  the  cell  length  and  repre¬ 
sent  the  highest  contributions  to  cell  voltage  losses.  Inversely,  the 
anode  concentration  overpotential  is  found  to  be  very  low.  This  can 
be  attributed  to  the  air  ratio  that  is  sufficiently  high  to  avoid  a  strong 
02  enrichment  of  the  anodic  gas. 

Fig.  9  gives  the  temperature  profiles  along  the  cell  length  calcu¬ 
lated  at  the  cathode/electrolyte  interface  and  in  the  interconnect 
plate.  In  the  SRU,  temperatures  are  found  to  vary  between  852  and 
865  °C.  This  temperature  range  is  consistent  with  the  exothermal 
operation  mode  chosen  in  this  reference  case  ( Ucen  =  1.5  V).  It  can  be 
noticed  that  the  temperature  profiles  are  approximately  symmet¬ 
rical  along  the  channel  length.  This  temperature  evolution  is  due 
to  the  heat  accumulated  into  the  cathode  and  anode  gases  flowing 
in  opposite  direction  (i.e.  in  a  counter-flow  configuration).  More¬ 
over,  the  heat  fluxes  removed  from  the  solid  structure  are  mainly 
exhausted  by  radiation  symmetrically  on  each  side  on  the  SRU. 
It  can  also  be  noticed  that  the  temperature  maximum  is  slightly 
shifted  from  the  middle  of  the  cell  towards  the  steam  inlet.  This 
phenomenon  is  consistent  with  the  extent  of  the  electrochemical 
reaction  along  the  cell  length:  the  reaction  is  favoured  at  the  steam 
inlet  where  the  current  density  is  the  highest  (cf.  Fig.  7). 


H2  outlet 
Air  inlet 


Fig.  7.  Evolution  of  the  local  current  density  along  the  cell  length  in  the  refer¬ 
ence  case  (Uceu  =  1.5  V,  steam  conversion  rate  =  96.3%,  condition  I  cell  gas  feeding, 
<5 electrolyte  =  90  |jim,  8anode  =  8cathode  =  50  |xm,  cell  active  area  =  77.44  cm2 ). 
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Fig.  9.  Temperature  evolution  plotted  as  a  function  of  the  cell  length  in  the  refer¬ 
ence  case  (Uceu  =  1.5  V,  steam  conversion  rate  =  96.3%,  condition  I  cell  gas  feeding, 
temperature  of  the  interconnect  plate  taken  on  the  surface  connecting  two  adjacent 
SRUs,  8 electrolyte  =  90  |jim,  8an0de  =  <5 cathode  =  50  fxm,  cell  active  area  =  77.44  cm2 ). 
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Fig.  10.  Polarisation  curves  obtained:  (i)  simulated  in  the  isothermal  state 
(T=800°C)  and  (ii)  simulated  with  the  connected  electrochemical  and  thermal 
model  taking  into  account  thermal  variations  with  operation  modes  (reference 
case,  condition  I  cell  gas  feeding,  8eiectroiyte  =  90  (Jim,  8anode  =  8cathode  =  50  (Jim,  cell  active 
area  =  77.44  cm2). 

The  temperature  of  the  interconnect  plate  is  found  to  be  close 
to  the  cell  temperature  (Fig.  9).  This  low  thermal  gradient  in  the 
stacking  direction  (y  direction)  is  due  to  the  adiabatic  assumption 
considered  here  for  the  modelling. 

3. 1.3.  Effect  of  cell  polarisation 

Fig.  10  provides  two  polarisation  curves  simulated  in  the  refer¬ 
ence  case.  These  curves  have  been  computed  either: 

(i)  by  considering  a  “virtual”  isothermal  operating  state  at 
T  =  800  °C  (i.e.  by  neglecting  any  thermal  effect), 

(ii)  or  by  coupling  the  thermal  and  electrochemical  models.  The 
thermal  variations  are  then  associated  with  endothermic  and 
exothermic  operation  modes.  In  this  case,  the  variation  of  cell 
temperature  with  the  current  density  (or  cell  voltage)  is  given 
in  Fig.  11. 


Fig.  11.  Cell  temperature  plotted  as  a  function  of  the  current  density  and  ther¬ 
mal  sources  (reference  case,  condition  I  cell  gas  feeding,  the  temperature  has  been 
taken  in  the  middle  of  the  cell  at  the  cathode/electrolyte  interface,  8eiectroiyte  =  90  |jim, 
8anode  =  8cathode  =  50  |jim,  cell  active  area  =  77.44  cm2 ). 


When  the  cell  operates  at  the  thermo  neutral  voltage  (-1.28  V), 
heat  sources  released  in  the  cell  such  as  ohmic  losses  are  counter¬ 
balanced  by  the  endothermic  nature  of  the  electrolysis  reaction.  As 
a  consequence,  the  cell  temperature  is  found  to  be  equal  to  800  °C 
(Fig.  1 1 ).  In  this  operating  condition,  the  current  density  is  found  to 
be  the  same  as  the  one  computed  in  the  “virtual”  isothermal  state 
(-0.68  A  cm-2  at  UceU  =  1 .28  V)  (Fig.  1 0). 

Below  the  thermo  neutral  voltage,  heat  sources  do  not  coun¬ 
terbalance  the  endothermic  electrolysis  reaction  and  the  cell 
temperature  is  found  to  be  lower  than  800  °C  (Fig.  11).  As  a  con¬ 
sequence,  the  cell  performance  calculated  in  this  case  is  lower  than 
the  one  obtained  in  the  isothermal  assumption.  Flowever,  the  dif¬ 
ferences  between  the  two  curves  remain  extremely  slight. 

Inversely,  above  the  thermo  neutral  voltage,  the  cell  warming 
and  its  impact  on  the  cell  polarisation  curve  are  more  pronounced. 
However,  it  has  been  calculated  that  the  activation  and  concen¬ 
tration  overpotentials  are  not  significantly  impacted  by  the  cell 
warming.  As  a  consequence,  the  cell  performance  improvement 
observed  above  the  thermo  neutral  voltage  cannot  be  explained  by 
the  activation  or  concentration  overpotentials  decrease.  Actually, 
this  phenomenon  must  be  attributed  to  the  decrease  of  the  elec¬ 
trolyte  ionic  resistance  with  the  temperature  increase.  This  effect 
is  illustrated  in  Fig.  1 2 :  above  the  thermo  neutral  voltage,  the  ohmic 
losses  calculated  in  the  “virtual”  isothermal  state  are  much  higher 
than  the  ones  simulated  by  coupling  the  thermal  and  electrochemi¬ 
cal  models.  It  can  be  seen  in  Fig.  12  that  the  difference  in  the  ohmic 
losses  between  the  two  cases  is  A(Rohmi)  «  40mV  at  1.0  A  cm-2. 
This  value  corresponds  roughly  to  the  cell  polarisation  curve 
improvement  calculated  at  the  same  current  density  in  Fig.  10. 

3.2.  Sensitivity  analysis:  effect  of  geometrical  and  operation 
parameters  on  steam  electrolysis 

3.2. 1.  Impact  of  hydrogen  over  steam  ratio  introduced  at  the 
cathode  side 

The  effect  of  cathodic  gas  composition  has  been  investi¬ 
gated  considering  the  Electrolyte  Supported  Cell  configuration 
anode  ~  8  cathode —  50  p>m  and  8  electrolyte =  90  |xmj.  The  SRU  perfor¬ 
mances  have  been  simulated  with  the  first  and  second  conditions 
of  gas  feedings,  for  which  the  H2  over  H20  ratio  R  is  decreased  from 
1/1  to  1/9  (see  Table  4). 

The  steam  molar  flow  rate  is  kept  constant  for  both  simulated 
cases.  That  condition  allows  obtaining  identical  SC  rates  accord¬ 
ing  to  the  applied  current  densities.  This  remark  means  that  the 
overpotentials  induced  by  the  steam  depletion  along  the  cell  are 
identical  in  both  simulations.  It  can  be  also  mentioned  that  the  total 
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Fig.  12.  Overpotential  plotted  as  a  function  of  current  density  and  SC  rate  calcu¬ 
lated:  (i)  in  the  isothermal  state  for  the  dashed  line  and  (ii)  by  coupling  the  thermal 
and  electrochemical  models  for  the  solid  curves  (reference  case,  condition  I  cell  gas 
feeding,  8eiectroiyte  =  90  pan,  8anode  =  8cathode  =  50  (Jim,  cell  active  area  =  77.44  cm2 ). 
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Fig.  13.  Polarisation  curves  simulated  for  conditions  I  and  II  of  gas  feedings 
(8 electrolyte  =  90  p.m,  8anode  =  8cathode  =  50  iJim,  cell  active  area  =  77.44  cm2 ). 


At  the  thermo  neutral  voltage  and  considering  the  same  operat¬ 
ing  conditions,  the  CSC  geometry  exhibits  better  performances  than 
the  ESC  configuration.  This  behaviour  is  obviously  induced  by  the 
difference  in  the  ionic  resistances  between  the  two  kinds  of  cells. 
However,  at  the  highest  investigated  cell  voltages,  it  can  be  seen 
that  the  gap  in  performances  between  CSC  and  ESC  tends  to  vanish. 
Indeed,  the  cell  temperature  elevation  decreases  the  YSZ  ionic  resis¬ 
tivity  and  its  relative  contribution  to  the  overall  cell  overpotentials. 
Moreover,  the  thick  cathode  substrate  of  the  CSC  geometry  hin¬ 
ders  the  diffusion  trough  the  electrode.  The  resulting  mass  transfer 
limitation  induces  large  concentration  overpotentials  and  leads  to 
the  limiting  current  density  observed  on  the  polarisation  curve 
(Fig.  15).  For  instance,  at  i  =  1 .2  A  cm-2  and  SC  =  90%,  the  term  p^ath 
is  increased  from  ~61  mV  for  ESC  to  ~228  mV  for  CSC  (considering 
the  first  condition  of  gas  feeding).  In  this  last  case,  the  drop  of  H2 
or  H20  molar  fraction  across  the  thick  cathode  becomes  significant 
as  shown  in  Fig.  16. 

3.2.3.  Impact  of  gas  dilution  through  the  cathode 

Considered  the  CSC  geometry,  the  high  concentration  overpo¬ 
tentials  have  been  correlated  to  a  mass  transfer  limitation.  This 


molar  flow  rate  and  air  ratio  are  maintained  to  the  same  levels  in 
such  way  that  the  thermal  conditions  (in  terms  of  heat  convec¬ 
tive  flow)  are  similar  in  both  simulated  cases.  The  balance  when 
decreasing  H2  molar  flow  rate  is  achieved  by  increasing  N2  molar 
flow  rate  (i.e.  by  increasing  the  dilution  term  §). 

A  significant  difference  has  been  found  between  the  two  simu¬ 
lated  polarisation  curves  (Fig.  13).  As  it  is  expected  by  considering 
Eq.  (12),  the  OCV  at  R  =  1/1  is  found  to  be  higher  than  the  one  cal¬ 
culated  at  R  =  1  / 9.  Whatever  the  investigated  current  density,  it  can 
be  seen  that  the  Area  Specific  Resistance  (ASR)  is  lower  for  R  =  1/1 
than  for  R  =  1  / 9.  The  contributions  to  the  voltage  losses  have  to  be 
analysed  to  understand  this  result: 

•  As  expected,  the  SRU  temperature  fields  have  been  found  to  be 
close  to  each  other.  As  a  consequence,  the  voltage  losses  due  to 
the  ohmic  resistances  are  quite  similar  in  both  simulated  cases. 

•  Since  the  condition  of  air  supply  is  the  same  for  both  polarisation 
curves,  the  anodic  activation  and  concentration  overpotentials 
are  found  to  be  not  impacted  by  the  modification  of  the  cathode 
gas  composition. 

•  Finally,  the  difference  between  the  two  polarisation  curves  is 
found  to  be  mainly  induced  by  the  increase  of  cathode  activa¬ 
tion  overpotential,  pcafth  (Fig.  14a),  and  in  a  lower  extent,  by  the 
increase  of  cathode  concentration  overpotential  p^ath  (Fig.  14b). 
The  impact  of  the  cathode  gas  composition  onto  is  ascribed 
to  the  dependence  of  the  exchange  current  density  i0  with  hydro¬ 
gen  and  steam  molar  fraction  (see  Eq.  (17a)).  Therefore,  it  can  be 
inferred  from  this  result  that  the  activation  overpotentials  can  be 
increased  from  negligible  values  to  substantial  ones  according  to 
the  condition  of  cathodic  gas  supply. 

3.2.2.  Impact  of  geometric  cell  configuration 

Electrolyser  performances  have  been  assessed  considering  a 
stack  assembled  with  Cathode  Supported  Cells  (CSC)  fabricated 
with  the  same  materials  as  taken  previously  in  the  ESC  case.  The 
simulations  have  been  carried  out  for  a  thickness  of  the  cermet 
substrate  equal  to  1000  p.m.  Thin  electrolyte  and  anode  layers 
deposited  onto  the  cathode  support  have  been  fixed  respectively  to 
15  p.m  and  50  p.m.  The  polarisations  curves  have  been  computed 
with  the  first  and  third  conditions  of  gas  feedings  (see  Table  4). 
These  simulated  curves  are  compared  in  Fig.  15  to  the  one  calcu¬ 
lated  in  the  reference  case  (i.e.  for  the  Electrolyte  Supported  Cell 
and  for  the  first  condition  of  gas  feeding). 
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Fig.  14.  Cathode  overpotentials  computed  with  conditions  I  and  II  of  gas  feedings 
(5 electrolyte  =  90  fxm,  8anode  =  8cathode  =  50  |jim,  cell  active  area  =  77.44  cm2 ):  (a)  activation 
overpotential  (b)  concentration  overpotential 
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Fig.  15.  Polarisation  curves  simulated  for  conditions  I  or  III  of  gas  feedings:  solid  line:  CSC  (8cathode  =  1000  (Jim,  8eiectroiyte  =  15  |jim,  8anode  =  50  (Jim,  cell  active  area  =  77.44  cm2). 
Dashed  line:  ESC  (8eiectroiyte  =  90  (Jim,  8anode  =  8cathode  =  50  (Jim,  cell  active  area  =  77.44 cm2). 


effect  is  found  to  be  worsened  when  the  steam  and  hydrogen  mix¬ 
ture  is  diluted  into  a  neutral  gas.  For  instance,  the  limiting  current 
density  is  roughly  decreased  from  ~1.2Acm-2  to  ~0.84Acnrr2, 
when  the  dilution  rate  £  is  raised  from  the  first  to  the  third  condi¬ 
tion  of  gas  feeding  (Fig.  15).  As  a  consequence,  addition  of  diluents 
into  H2,  H20  mixture  leads  to  a  decrease  in  the  cell  performances. 

This  effect  of  diluent  concentration  on  cell  limiting  current  is 
consistent  with  the  results  obtained  in  SOFC  mode  by  Jiang  and 
Virkar  [22].  These  authors  have  shown  that  the  maximum  achiev¬ 
able  current  density  by  the  fuel  cell  is  limited  at  high  concentrations 
of  diluent.  This  phenomenon  is  explained  by  the  hydrogen  and 
steam  effective  gas  diffusivity  which  is  lowered  with  increasing  the 
nitrogen  concentration.  Indeed,  it  can  be  easily  demonstrated  with 


Fig.  16.  Evolution  of  the  H2  and  H20  molar  fraction  along  the  CSC  cell 
length  (Uceu  =  1.5  V,  steam  conversion  rate  =  90%,  condition  I  of  cell  gas  feeding, 
8Cathode  =  1 000  |Jim,  8eiectroiyte  =  15  |jim,  <5 anode  =  50  (Jim,  cell  active  area  =  77.44  cm2 ). 


Eqs.  (Al-2)  and  (Al-3)  that  the  FI2  and  H20  effective  diffusivities 
depend  on  the  N2  concentration  as  followed: 
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It  can  be  noticed  that  the  FI20/N2  binary  effective  diffusion  coeffi¬ 
cient  is  lower  than  the  one  for  H20/H2  species  (D^0  n2  =  2.858  x 

1CT5  m2  s-1  <  DeJfn  H  =  1.004  x  10-4  m2  s-1).  In  this  case, 

h2o,h2  j 

according  to  Eq.  (27),  the  term  D^0  is  decreased  with  increasing 
the  concentration  of  diluent  yN2.  This  trend  is  illustrated  in  Fig.  17 
where  the  steam  effective  diffusion  coefficient  has  been  plotted  as 
function  of  nitrogen  molar  fraction. 


4.  Discussion 


4.1.  Effect  of  radiative  heat  losses  on  the  stack  temperature 

As  shown  in  Fig.  11,  the  SRU  maximum  temperature  reaches 
roughly  865  °C  at  1 500  mV  and  for  i  ~  1 .28  A  cm-2.  In  this  case,  it  is 
worth  underlining  that  the  main  part  of  heat  sources  is  removed 
by  the  radiation  losses.  Indeed,  73%  of  the  heat  generated  by  the 
electrochemical  process  within  the  cell  is  dissipated  by  radiative 
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Fig.  17.  Effective  steam  diffusion  coefficient  D^Q  plotted  as  a  function  of  N2  molar 
fraction  (7=800°C,  £  =  0.4,  r  =  1  /xm,  PT=  1  atm). 


transfer  between  the  SRU  and  the  insulating  envelope.  Only  27% 
of  the  remaining  heat  sources  are  transferred  to  anode  and  cath¬ 
ode  gas  streams.  It  can  be  inferred  from  this  result  that  the  radiative 
losses  have  a  determinant  impact  on  the  cell  temperature:  the  radi¬ 
ation  between  the  stack  and  its  insulating  envelope  appears  to  be 
the  most  efficient  thermal  transfer  mode  to  evacuate  (or  bring) 
the  heat  generated  (or  absorbed)  by  the  electrolysis  process.  Con¬ 
versely,  the  convective  heat  transfer  remains  limited  along  the  gas 
channels.  As  a  consequence,  the  radiation  represents  the  main  heat 
transfer  mode  and  it  cannot  be  neglected  in  the  steam  electrolysis 
modelling. 

Because  of  its  significant  importance,  the  stack  temperature 
should  strongly  depend  on  the  radiation  efficiency.  To  assess  this 
effect,  a  sensitivity  analysis  has  been  carried  out  by  varying  the 
factor  A  which  links  the  radiative  losses  to  the  solid  temperatures: 

A  =  f.  d<p2/  (t-4  -  r4nsutol.ng)  (see  Eq.  (25)).  In  Fig.  18,  the  simulated 

cell  temperature  is  plotted  as  a  function  of  this  factor  A.  The  percent¬ 
ages  of  radiative  and  convective  heat  losses  have  been  also  added 
into  the  figure.  It  can  be  seen  that  the  relative  proportions  of  heat 
removed  by  radiation  and  convection  remain  in  the  same  order  of 
magnitude  whatever  the  value  of  A  (since  the  thermal  sources  are 
equal  for  each  simulated  points).  When  the  “efficiency”  of  the  radi- 
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Fig.  18.  Cell  temperature  plotted  as  a  function  of  A  (Uceu  =  1.5  V,  Condition  I  cell 
gas  feeding,  the  temperature  has  been  taken  in  the  middle  of  the  cell  at  the 
cathode/electrolyte  interface,  Seiectroiyte  =  90  p,m,  8anode  =  8cathode  =  50  p,m,  cell  active 
area  =  77.44  cm2). 


ation  is  lowered  by  decreasing  A,  the  cell  temperature  is  observed  to 
be  strongly  increased.  Indeed,  the  solid  temperature  has  to  be  raised 
to  maintain  a  sufficient  radiative  heat  flow  and,  hence,  to  evacuate 
the  heat  generated  within  the  cell.  This  phenomenon  maintains 
the  energy  balance  between  the  heat  sources  and  the  radiation  and 
convection  fluxes.  As  a  consequence,  a  limitation  of  the  radiation 
heat  exchange  into  the  electrolyser  stack  will  lead  to  an  increase 
of  the  SRU  temperature.  This  temperature  elevation  could  be  detri¬ 
mental  regarding  to  the  metal  interconnect  plates  or  the  contact 
layers.  Moreover,  the  SRU  ageing  processes  will  be  accelerated  at 
this  temperature  level.  In  this  case,  an  operating  condition  close  to 
the  thermo  neutral  voltage  will  be  the  only  way  to  keep  a  reasonable 
temperature  within  the  SRU  solid  structure. 

4.2.  Effect  of  activation  and  concentration  overpotentials  on  the 
electrolyser  performances 

4.2 A.  Activation  overpotential 

As  mentioned  in  Section  3.1 ,  a  substantial  contribution  to  the  cell 
polarisation  has  been  attributed  to  the  anode  activation  overpoten¬ 
tial  p™°de  (Figs.  8  and  12),  while  the  cathode  activation  polarisation 
pcafth  represents  only  a  small  part  of  the  cell  voltage  losses.  Although 
the  term  p^f1  can  be  significantly  affected  by  the  H2/H20  ratio  as 
shown  in  Section  3.1  (Fig.  14),  the  global  electrochemical  behaviour 
of  the  electrolyser  cannot  be  strongly  deteriorated  by  this  over¬ 
potential  (Fig.  13).  That  means  that  H2/H20  ratio  can  be  chosen 
according  to  other  criteria  such  as  cell  materials,  durability  or  sys¬ 
tem  efficiency. 

4.2.2.  Concentration  overpotential 

In  the  present  simulated  conditions,  the  anodic  concentration 
overpotential  Pconc6  contributes  less  than  20  mV  to  the  irreversible 
losses  (Fig.  12).  This  result  means  that  this  contribution  remains 
insignificant  whatever  the  investigated  operating  conditions  and 
cell  configurations. 

For  the  ESC,  the  cathodic  concentration  overpotential  pccf£c 
reaches  only  moderate  values  even  at  high  current  densities  and 
SC  rate.  For  instance,  this  overpotential  represents  approximately 
~1 1%  of  the  total  voltage  losses  at  i  =  1 .2  A  cm-2  and  SC  =  90%  (in  the 
reference  case). 

Conversely,  considering  the  CSC  design,  the  cathodic  concen¬ 
tration  overpotential  pccfdc  is  liable  to  provide  a  large  contribution 
to  the  cell  polarisation.  For  instance,  this  overpotential  represents 
up  to  ~41%  of  the  irreversible  losses  at  i=  1.2  A  cm-2  and  SC  =  90% 
(condition  I  of  gas  feeding  and  hcathode  =  1000  [xm).  Furthermore, 
the  Pcfnc  polarisations  are  exacerbated  when  diluents  are  added 
into  the  steam  and  hydrogen  mixture  (see  Section  3.2.3).  These 
high  cathode  concentration  overpotentials  have  been  ascribed  to 
a  mass  transfer  limitation  through  the  thick  cathode.  Because  of 
this  limitation,  the  difference  between  the  CSC  and  ECS  designs 
tends  to  vanish  at  high  current  density  and  SC  rate  (although  the 
CSC  configuration  exhibits  better  performances  than  ESC  geometry 
at  intermediate  cell  polarisation:  Fig.  15).  It  is  worth  underlining 
that  this  effect  will  depend  strongly  on  the  thickness  of  the  cath¬ 
ode  substrate.  As  it  is  illustrated  in  Fig.  19,  the  cell  limiting  current 
is  increased  and  its  associated  cathode  concentration  overpoten¬ 
tials  are  decreased  when  the  thickness  of  the  support  is  reduced. 
For  instance,  when  the  cathode  thickness  is  decreased  to  600  p,m, 
the  term  p^fjf  is  lowered  below  50  mV  at  0.8  A  cm-2  (Fig.  19b).  It 
can  be  noticed  that  this  last  value  is  consistent  with  the  simulation 
result  established  by  Udagawa  et  al.  [13]  in  similar  conditions. 

As  a  consequence,  the  thinning  of  the  cermet  substrate  improves 
the  cell  performances  but  also  reduced  its  mechanical  robustness. 
Therefore,  a  compromise  remains  to  be  found  from  this  point  of 
view.  Improvements  of  gas  transport  in  the  porous  cathode  could 
help  to  obtain  these  requirements. 
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Steam  conversion  rate  (%) 


Current  density  (A  cnv2) 

Fig.  19.  Effect  of  cathode  thickness  on  CSC  electrochemical  behaviour  (condition 
I  of  gas  feeding,  Seiectroiyte  =  15  |jim,  8anode  =  50  |jim,  cell  active  area  =  77.44 cm2),  (a) 
Polarisation  curves,  (b)  Cathodic  concentration  overpotential  plotted  as  a  function 
of  current  density. 


temperature  increase  of  the  SRU,  an  operating  condition  close  to 
the  thermo  neutral  voltage  is  necessary. 

Parametric  study  has  been  carried  out  to  analyse  the  SOEC 
irreversible  losses.  It  has  been  found  that  anode  activation  over¬ 
potentials  are  significant.  Conversely,  the  polarisation  due  to  the 
cathode  activation  is  much  more  restricted.  Thus,  even  if  this  acti¬ 
vation  overpotential  depends  on  the  steam  and  hydrogen  molar 
fraction  introduced  at  the  cathode  inlet,  the  H2/H20  ratio  can  be 
fixed  according  to  other  criteria  such  as  cell  materials,  durability  or 
system  efficiency. 

The  anode  concentration  overpotentials  have  been  found  to  be 
insignificant  whatever  the  operating  condition.  The  cathode  con¬ 
centration  overpotentials  are  found  to  be  moderate  for  ESC  design, 
while  they  can  increase  drastically  with  current  density  for  CSC 
geometry.  Finally,  addition  of  diluent  into  the  H2,  H20  mixture 
increases  these  overpotentials  since  the  effective  gas  diffusion  of 
steam  is  substantially  lowered. 
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Appendix  A.  Steam,  hydrogen  and  nitrogen  diffusion 
through  the  porous  cathode 

The  DGM  model  (see  Eq.  (3))  is  applied  for  a  gas  mixture  com¬ 
posed  of  hydrogen,  steam  and  nitrogen.  It  leads  the  three  following 
expressions: 


In  addition,  it  has  been  found  that  the  contribution  to  the  gas 
transport  of  both  molecular  and  Knudsen  diffusions  is  roughly  in 
the  same  order  of  magnitude.  Then,  if  one  of  these  diffusion  pro¬ 
cesses  is  neglected  in  the  modelling,  the  simulation  results  will 
underestimate  the  cathode  concentration  overpotential. 

5.  Conclusion 


(Al-2) 


A  2D  multi-physic  in-house-model  has  been  developed  to  anal¬ 
yse  the  performances  of  SOEC  stack.  The  model  encompasses  a 
combined  electrochemical  and  thermal  description  of  the  elec- 
trolyser.  In  order  to  facilitate  the  numerical  implementation,  an 
analytical  solution  for  multi-species  diffusion  across  the  porous 
cathode  was  derived  from  the  DGM  equations. 

Understandable  predictions  of  SOEC  behaviour  can  be  achieved 
by  using  the  model  since  it  provides  all  the  important  parameters  of 
the  stack  operation:  distribution  of  local  temperature,  heat  fluxes, 
current  densities,  gas  concentrations  and  overpotentials. 

Simulations  have  shown  that  thermal  equilibrium  of  the  stack 
is  strongly  dependant  on  radiative  heat  losses  with  the  insulating 
envelope.  In  the  exothermic  operation  mode,  the  cell  temperature 
depends  on  the  radiation  efficiency:  in  order  to  avoid  a  too  high 


Eqs.  (Al-l)-(Al-3)  constitute  a  set  of  3  equations  including  3 
unknowns  (i.e.  yn2,yw2  and  yH2o)-  Since  neither  consumption  nor 
production  of  nitrogen  occurs  at  the  cathode  side  (meaning  that 
NN2  =  0),  Eq.  (Al-1),  can  be  written  as  followed: 

d^L=k,xyN2  (Al  -4) 
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Appendix  B.  Temperature  of  fluid  into  the  gas  manifold 
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This  last  differential  equation  can  be  solved  considering  the 
boundary  condition,  y^2  =  y™"al  when  y  =  8C: 

yw2  =  yCN2nal  x  exp  {fc!  (y  -  Sc) }  (Al-5) 

The  expression  of  the  molar  fraction  of  nitrogen  at  the  cath¬ 
ode/electrolyte  interface  is  directly  deduced  from  (Al-5): 

=  y™2al  x  exP  (_ki Sc)  (A1 

Assuming  that  no  accumulation  of  species  occurs  in  operation 
(i-e.  Nh2  +Nh2o  =  0)  and  withyH2  +yH2o  +yn2  =  1.  Eq.  (Al-2)  can 
be  written: 

k2  +  k3xyN2  =  -d^  (Al-7) 

with 


fo  = 


k  3  = 


RxT  /  Nh2 


Nh2 


PT  \ 

Deff 

\UH2,k 

Deff 

^h2,h2o 

RxT  | 

(  Nh2 

Nh2 

PT  \  Deff 
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0 


h2,n2  h2,h2o 


The  boundary  condition  of  Eq.  (Al  -7)  corresponds  toyn2  =  y™™1 
wheny  =  <5c.  The  solution  of  this  differential  equation  corresponds 
to: 


yH2  =  yf.7'  +  Me  +  ^yi7'  -  M -  |ycNQ2na' exp  {My- sc)} 

(Al-8) 


The  aim  of  this  section  is  to  determine  the  gas  temperature  evo¬ 
lution  along  the  manifold  at  the  anode  or  cathode  inlets  (Fig.  1). 
These  gas  temperature  evolutions  Tg(z)  are  determined  by  writing 
an  energy  balance  between  z  and  z  +  dz: 

n(z)Cp{Tg(z  +  dz)-Tg(z)}  =2(a  +  b)dzxhx  {Ts-Tg(z)}  (A2-1) 

where  Cp  and  n(z)  are  respectively  the  specific  heat  and  the  molar 
flow  rate  of  the  gas  mixture.  The  terms  a  and  b  denote  respectively 
the  manifold  width  and  height.  As  illustrated  in  Fig.  3,  the  molar 
flow  rate  n(z)  depends  on  the  z  coordinate  along  the  canal: 

n(z)  =  ntotx{  1-5}  (A2-2) 

where  ntot  is  the  total  molar  flow  rate  introduced  into  the  SRU  and 
L  the  manifold  length. 

Eqs.  (A2-1 )  and  (A2-2)  lead  to  the  following  first-order  differen¬ 
tial  equation: 

ntot  { 1  “  I }  Cp  x  +  h2  (a  +  b)  X  Tg(z )  =  h2(a  +  b)Ts  (A2-3) 

This  last  equation  can  be  solved  assuming  that  the  wall  temper¬ 
ature  of  the  solid  phase  Ts  is  constant  in  the  z  direction  (cf.  model 
assumption  in  Section  2.1): 

r  r  Z  1  1  ((h2  (a+b)L)/ntotCp) 

Tg(z)  -Ts  ex  yitot  |l  -  ijj  (A2-4) 

At  the  SRU  inlet,  Tg(z  =  0)  =  Tintet  =  800 0  C : 

Tg(z)  =  (Tinlet  -  Ts)  x  [nt0 tf-h2(a+b)L)/ntotCp) 
r  r  Z  1  1  ((h2  (a+b)L)/ntotCp) 

x{ntotY~L}\  +Ts  (A2_5) 
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At  the  cathode/electrolyte  interface,  this  expression  becomes: 


y|?2  =  yCH™1  +  k2sc  +yLNa2nal |  { l  -  exp (-k, 5C)} 
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Considering  the  same  previous  assumptions,  Eq.  (Al-3)  can  be 
also  written: 
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Integration  of  this  last  differential  equation  allows  determining 
the  steam  molar  fraction  at  the  cathode/electrolyte  interface: 

y|?2o  =yH2o'  +  Mc+y|71||  {i  -exP(-k,5c)}  (ai-ii) 

It  is  worth  noting  that  the  constants  /<i ,  /<2,  k3,  k4  and  can  be 
expressed  as  function  of  the  current  density  by  using  Eq.  (9).  There¬ 
fore,  Eqs  ( Al  -6),  (Al  -9)  and  ( Al  -1 1 )  can  be  used  to  calculate  the  gas 
composition  at  the  cathode/electrolyte  interface  as  a  function  of  the 
molar  fractions  in  the  channel  and  the  local  current  density. 
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